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An improved hybrid particle-finite element method has been developed for the simulation of 
hypervelocity impact problems. Unlike alternative methods, the revised formulation computes 
the density without reference to any kernel or interpolation functions, for either the density 
or the rate of dilatation. This simplifies the state space model and leads to a significant 
reduction in computational cost. The improved method introduces internal energy variables 
as generalized coordinates in a new formulation of the thermomechanical Lagrange equations. 
Example problems show good agreement with exact solutions in one dimension and good 
agreement with experimental data in a three dimensional simulation. 
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INTRODUCTION 

Studies of hypervelocity impact phenomena are motivated by a variety of science and en- 
gineering applications [1]. Examples include scientific research on planetary impacts [2] and 
equations of state [3] and engineering research on the design of spacecraft shielding [4] and 
kinetic energy penetrators [5] . The proceedings of a recent international symposium [1] show 
that the use of computer simulation in this field is increasing, as improvements in numerical 
methods and computing power make it possible to address problems of greater complexity 
and larger scale. Simulation is of particular importance, as an adjunct to experimental work, 
when material costs are high [6] or when impact velocities beyond the range of light gas guns 
are of interest [7]. 

Simulation work in this field has applied a number of different numerical methods, based 
on continuum mechanics, particle dynamics, or mixed kinematic schemes. Continuum meth- 
ods [8] employ either an Eulerian hydrodynamic [9,10] or a Lagrangian finite element [11] 
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approach, or some Arbitrary Lagrangian-Eulerian (ALE) based generalization of these tech- 
niques [12,13]. A large majority of particle codes employ a smooth particle hydrodynamics 
(SPH) technique [14,15,16], although some alternative particle based methods have been 
proposed [17]. Some disadvantages of pure continuum or pure particle based methods [18] 
have motivated the development of mixed continuum-particle formulations [4,19,20]. The 
most widely used mixed method is a coupled particle-finite element technique [19]. This 
technique initializes distinct material regions with either SPH particles or or Lagrangian 
finite elements, then quantifies subsequent material interactions using a particle-to-surface 
contact-impact algorithm. The alternative coupled particle-element method of Johnson and 
co-workers [20] maps damaged or failed elements into particles, and again quantifies particle- 
element interaction using a special contact algorithm. Both these methods are subject to 
tensile instability and numerical fracture problems. 

The alternative mixed method of Fahrenthold and co-workers [21] is based on a hybrid 
particle-finite element formulation. This method in not subject to tensile instability and 
numerical fracture problems and eliminates the requirement for special treatment of particle- 
to-element contact-impact. It avoids both the mass diffusion problems of Eulerian methods 
and the mass and energy discard associated with Lagrangian element erosion algorithms. It is 
labeled a hybrid (versus coupled) method since it introduces both elements and particles for 
all material control volumes, then employs the elements and particles in tandem to represent 
distinct physics. The particles model all inertia, contact-impact, and thermomechanical 
response in compressed states, while the elements model tension and elastic-plastic shear. 
The method incorporates both ellipsoidal particles and time varying particle volumes and as a 
result can represent large density variations with relatively small neighbor counts. Previous 
work employing nonspherical evolving kernels has been rather limited and most particle 
simulations represent high densities using spherical particles, a fixed contact length, and 
relatively large neighbor sets. 

The preceding formulation combines a true Lagrangian description of material strength 
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effects with a general particle based model of contact-impact dynamics, and has been vali- 
dated in simulations of impact experiments conducted at velocities ranging from one to ten 
kilometers per second [22], In the hypervelocity impact regime, where large strain plastic- 
ity, perforation, fragmentation, melting, and complex multi-structure contact-impact effects 
are often present, this formulation provides a particular combination of advantageous fea- 
tures not offered by alternative numerical methods. The present paper describes an im- 
proved hybrid particle-finite element method, modifying the formulation of Ravishankar 
and Fahrenthold [21] in two respects. First, the density is determined by integrating non- 
holonomic constraints imposed on the system level thermomechanical model. Second, the 
entropy states used previously to model the thermal domain are replaced by particle inter- 
nal energies. These modifications simplify the method, reduce its computational cost, and 
incorporate equations of state expressed in standard functional form. 

Unlike alternative methods, the revised formulation eliminates entirely the use of kernel 
or interpolation functions to represent the density or rate of dilatation fields. The density 
evolution equations are developed by direct reference to large deformation kinematics, avoid- 
ing any requirement to specify the functional dependence of an interaction potential on the 
particle coordinates. The latter task has proven to be quite difficult in an SPH context and is 
a principal focus of the general particle dynamics literature. The revised method introduces 
the use of internal energy variables as generalized coordinates in a new thermomechanical 
formulation of the discrete Lagrange equations. This avoids the requirement to construct 
Legendre transforms of the internal energy function, in order to express the dependence of 
pressure and temperature on entropy. 

The present paper is organized as follows. First the particle and element kinematics are 
defined, followed by the kinetic co-energy and thermomechanical potential energy functions 
for the particle-element system. Second the evolution equations for the density are developed, 
followed by the evolution equations for the plastic and damage variables, all of these relations 
representing nonholonomic constraints on the system level model. Third the numerical 
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viscosity and numerical heat diffusion models are introduced, and evolution equations for 
the internal energy state variables are described, the latter states serving as generalized 
coordinates in a thermomechanical Lagrangian formulation. Fourth the discrete Lagrange 
equations for the particle-element system are derived, taking an explicit state space form 
convenient for numerical implementation. Finally application of the method is illustrated 
in one dimensional problems with exact solutions and in three dimensional simulations of 
representative hypervelocity impact problems. 


PARTICLE KINEMATICS 

The inertia of the modeled system is represented by a collection of n ellipsoidal particles, 
with the mass of particle i and h^\ h^\ the half-lengths of its major axes. The 
position and orientation of each particle is determined by its center of mass position vector 
( cW ) and an Euler parameter [23,24] vector ( ) 
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where a superscript T denotes the transpose. 

It is convenient to note here certain properties of Euler parameters, and to cite a number 
of well known [23] kinematic relations associated with their use. The Euler parameters 
provide a singularity free description of arbitrary particle rotations. They define a rotation 
matrix (R^) for each particle 
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which relates vector components v described in a fixed global Cartesian coordinate system 
to corresponding components v described in a co-rotating system aligned with the particle 
major axes, using 

v = R w v (5) 

The Euler parameters and their time derivatives are related to the angular velocity vector 
of the particle ( c*/ 1 ' ), described in the co-rotating frame, by 

e w = ^ G wr u> (i) (6) 

Similarly the antisymmetric matrix with axial vector which satisfies 

O w v = w (i) x v (7) 

for all vectors v, is related to the Euler parameters and their time derivatives by the relations 

n (i) = 2 G (i) G (i)T = -2 G (i) G (i)r = R (l)T R (i) (8) 


For ellipsoidal particles it is convenient to describe the separation distance of the mass 
centers for particles i and j using the ellipsoidal coordinate 


C (iJ) = (c (i) - c 


0)) T Hd) ( c (0 _ c^) 


(9) 


defined in the co-rotating system of particle j using 
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where the constant [3 allows for close packing at the reference density. The time derivative 
of this ellipsoidal coordinate , defined for i ^ j, is 
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where 


r (*j) _ c (*) _ c 0) 




(13) 


and may be used to quantify the rate of compression for an array of ellipsoidal particles. 

The preceding results will be used in later sections to account for rotational inertia and 
particle kinematics not present in the vast majority of particle models, which assume a 
spherical particle geometry. 

FINITE ELEMENT KINEMATICS 

This section describes the finite element kinematics employed in the present paper. The 
elements used here are eight noded hexahedra, well known and described in detail by Hal- 
lquist [11] and others. Since all inertia effects are represented by the particles, no mass 
matrix is defined. 

Each structure in the model is subdivided into uniform hexahedra with orthogonal faces, 
with ellipsoidal particles located at each node and at the centroid of each element. The 
center of mass coordinates for particles located at element vertices are also nodal coordi- 
nates for the hexahedra, and are used to compute the shearing strain. The center of mass 
coordinates for particles located at the element centroids are used, in combination with the 
nodal coordinates, to define six subelements for each hexahedron. The volumes of these 
subelements are used to compute interparticle tension forces. 

The following Lagrangian finite strain deformation measures [25] are used in the stored 
energy functions for the elements, associated with tension and shear, and in the plastic 
constitutive relations. The shear strain for element j is 



(14) 


where 



and F-T is an element deformation gradient computed using one point integration [11], The 
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clastic shear strain is defined as 


E e 0') — — E p(j ') 


(16) 


where (E p -^) is a plastic stain tensor whose flow rule satisfies the isochoric plastic deformation 
constraint 


tr = 0, K" ,J I = i - 1) (17) 


( 17 ) 


The subelement Jacobians are denoted by J^ ,k \ where the index k designates one of six 


subelements for the jth hexahedron. 

KINETIC CO-ENERGY AND POTENTIAL ENERGY 

An energy method (Lagrange’s equations) is adopted here, to facilitate the systematic 
integration of diverse particle and element based modeling concepts. The stored energy 
functions considered here are a kinetic co-energy function for the particles, an internal energy 
function for the particles, and element potential energy functions which account for tension 
and shear. Damage variables are introduced to model element failure in a thermodynamically 
consistent fashion. Constitutive assumptions different from those adopted here may be 
introduced without change to the underlying methodology. 

The system kinetic co-energy is the sum of the particle co-energies 


n 



(18) 


where T*W is the co-energy for particle i, due to translation and rotation 



2 2 
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with Jh) a constant moment of inertia matrix described in the co-rotating particle frame. 
The system kinetic co-energy function defines the generalized momenta 




( 20 ) 


where pW and are translational and angular momentum vectors for the Ah particle. 
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The system potential energy has the general form 


1/ = u w + y v o U) ^ U) + Y v o (j,k) ^ {j,k) 
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( 21 ) 


where is the total internal energy for particle i and the pressure (PW) and temperature 
($W) are described by an equation of state [26] with functional form 

pb) = p(0 (pW , U W), w (p w , ) 


( 22 ) 


with pW and u W the density and the internal energy per unit mass 
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The second term depends on the number of elements (n e ), the reference volume ( Vo ^) for 
element j, and the strain energy per unit volume in shear (V^), here assumed to be 


= (1 - d U) ) /i (i) tr (E e ^ T E e ^ 
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where d (j> is a shear damage variable and is a shear modulus. The third term depends 
on the number of subelements per element (n s ), the subelement reference volumes (Vo’®), 
and the strain energy per unit volume in tension (' 0 b») ; here assumed to be 


= - (1 - D®) K& < - 1 > 2 


(25) 


where D <3> is a normal damage variable, K (j> is a bulk modulus, < x > denotes the bracket 
function 

< x > = x u{x) (26) 

and u denotes the unit step function. Since the subelement Jacobians and the shear strain 
tensor depend on the particle center of mass coordinates 


= j(j,k) ( c (0) ) E e b) = E e(j) (c w , E ?j(j) ) 


it follows that the system potential energy has the general functional form 
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The system potential energy defines the generalized conservative forces 
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as well as the deviatoric stress 
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due to damage evolution. Note that when the internal energy is introduced as a generalized 
coordinate, the associated generalized forces are constant. 

With the system Lagrangian now defined, the next four sections describe evolution equa- 
tions for the internal state variables. 


DENSITY EVOLUTION RELATIONS 


This section derives density evolution relations for the particles, by extending certain 
exact results for the deformation kinematics of a unit cell of spherical particles, arranged 
in a body centered cubic packing scheme. For uniform compression of such a unit cell, in 
isolation, the cell center density pW is related to the reference density p 0 ^\ the center of 
mass separation distances and particle radii h (j) for its eight neighbors by 


pf 8 h 
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If the particles are not spherical but ellipsoidal then 
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This expression for an isolated cell may be extended to a cell array of arbitrary size by 
describing the same kinematics in rate form. Taking the time derivative of the last equation 
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and multiplying by a step function which allows for contact with near neighbors only 


yields 
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where the summation is now over all n particles. The coefficient bbbh) must satisfy 
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in order to correctly reflect the dependence of particle contact distance on the local density. 
Otherwise particles will contact remote neighbors through intervening matter. Hence 
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where denotes the Kronecker delta. Note that W^ 1 ’ 31 performs no interpolation. Intro- 
ducing the kinematic relation for Q l,3 \ developed in an earlier section, yields 
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which is the constraint form of the density evolution relations. The coefficients of the particle 
translational velocities and angular velocities in this expression will determine generalized 
forces in the momentum balance (Lagrange) equations derived in a later section. When the 
particle velocities in this expression are eliminated in favor of the particle momenta 


c « = m (i)- 1 p (i) , 




( 37 ) 


the density evolution equations take an explicit state space form convenient for use in nu- 
merical simulation. 


PLASTICITY AND DAMAGE MODELS 

This section introduces evolution equations for the plastic and damage variables. As in the 
case of the potential energy, alternative constitutive assumptions may be introduced without 
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change to the basic modeling methodology. The plastic flow rule used here is adapted from 
reference [27], and represents the simplest possible accommodation of the aforementioned 
isochoric plastic deformation constraint. The flow rule is 
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and the invariant operator is defined by 
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for any second order tensor T. The fourth order tensor coefficients in the flow rule are 
defined by 
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for any symmetric second order tensor T. The yield function is 
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where Y (j > is the yield stress 
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with e p b) the effective plastic strain, k ( j> a strain hardening coefficient, a b) a strain hardening 
exponent, 77 b) a thermal softening coefficient, and 9 H b) the homologous temperature. The 
effective plastic strain is determined by integrating the rate relation 
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while the incremental plastic strain for a time step At is computed using 


*><,-) <||s»«||-y«> 
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The damage evolution equations applied here are adapted from reference [28], and dissi- 
pate the strain energy stored in tension and shear over fi time steps, once an element meets 
any stipulated material failure criteria. The evolution equations are 


D (*> = u( 1 - D U) ), 
hAt K h 
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where is initialized to zero, and is set to a value of one when the accumulated plastic 
strain, temperature, or element compression reach corresponding critical values for the plastic 
failure strain (e^ ), melt temperature (dm), or maximum compression (JrP)- Other failure 
criteria may of course be specified. 

In general terms, the plastic and damage evolution equations are noholonomic constraints 
of the form 
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on the system level Lagrangian model. 


ARTIFICIAL VISCOSITY AND HEAT DIFFUSION 

Shock physics codes of the continuum or particle type incorporate a numerical viscosity 
and artificial heat diffusion. The forms used here are typical of particle codes, with one 
exception. Since the ellipsoidal particles used here admit rotational degrees of freedom, a 
viscous torque has been added which damps the relative rotation of neighboring particles. 

A viscous force is introduced for converging particles only 
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where the relative normal velocity is 
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and the viscosity coefficient is 
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with Cs'^ and V} 1 ' 1 a sonndspeed and particle reference volume. The parameters c Q and c\ are 
nondimensional linear and quadratic numerical viscosity coefficients. 

Similarly the viscous torque is 
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Finally the thermal power flow due to artificial heat diffusion is taken to be 
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where the heat transfer coefficient is 
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with ct a specific heat and k 0 a numerical heat diffusion coefficient. 


INTERNAL ENERGY EVOLUTION EQUATIONS 

The last internal state variable to be considered is the internal energy. The introduction 
of internal energy states as generalized coordinates allows the thermomechanical problem of 
interest here to be solved using energy methods. 

The internal energy evolution equations for particle i are 

lj(i) — jjwrk(i ) _|_ jjirr(i ) _ jjcon{i ) Pggh 
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where the first term accounts for mechanical work, the second term accounts for the effects 
of irreversible entropy production, and the third term represents numerical heat diffusion. 
The mechanical power flow for particle i is 
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The energy dissipation due to irreversible entropy production for particle i depends on the 
viscous forces and torques, which act on the particles, and on the dissipation in the elements 
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where Q irr(j > is a power flow due to damage evolution and plastic deformation in element j 
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and is the fraction of the dissipation in element j associated with particle i. 

Finally the internal energy flows due to numerical heat diffusion are 
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As in the case of the density evolution equations, a constraint form of the internal evolution 
relations is used to identify the generalized forces which appear in the Lagrange equations 
developed in the next section. For numerical implementation of the method, the generalized 
velocities are eliminated by introducing the momentum states as well as evolution relations 
for the density, plastic, and damage state variables. The resulting internal energy evolution 
relations take an explicit state space form. 


LAGRANGE’S EQUATIONS 

The preceding sections defined stored energy functions and nonholonomic constraints for 
the thermomechanical particle-element system. This section develops the final ODE model. 
The results of Shivarama and Fahrenthold [21] allow in the present case Lagrange’s equations 


to take the canonical form 
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where q c W, Q C W, Q p W, Q u< ^\ Q d ^\ Q D< ^\ and Q p V) are generalized forces determined by 
the nonholonomic constraints. The degenerate forms of the Lagrange equations for the 
internal state variables are due to the fact that those variables are not associated with any 
generalized momenta. Introducing Lagrange multipliers 7 p h) ; ^u(i) ^ y(j) , ^D{j) ^ arK ] 
for the constraints, the generalized forces are found to be 
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These results allow the unknown Lagrange multipliers to be determined in closed form, so 
that the final Lagrange equations are 
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where the generalized forces and torques due to particle interactions are 
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Supplemented by the evolution equations for density, internal energy, shear damage, normal 
damage, and plastic strain, the result is an explicit first order ODE model for the thermo- 
mechanical particle-element system. 


EXAMPLE SIMULATIONS 

This section describes four example problems which illustrate application of the improved 
particle-element method developed in this paper. The first two examples involve one dimen- 
sional test problems with known exact solutions. The third and fourth example problems 
involve three dimensional simulations. The third example models a published experiment 
while the fourth example measures the relative computational cost of the present method in 
a hypervelocity impact application of current research interest. 

The first example is the wall shock problem of Noh [31]. The simulations employ an ideal 
gas equation of state 

0 = — (u - u 0 ) 
c v 


P = (7 - 1) P (u - O, 


(80) 


and the parameters shown in Table 1. This problem models the collision of a fluid stream 
located initially in the region 0.0 < x < 0.5 with a rigid wall located at x — 0. The initial 
conditions are p = p D , u = u 0 , and v = —1. Figures 1 through 4 plot the simulation results for 
velocity, density, pressure, and temperature at the stop time of 0.3, for a model composed of 
200 particles. The numerical results show good agreement with the exact solution, although 
better results have been obtained using finite difference and finite element methods [31,35]. 
Table 2 shows convergence of the simulation results, as the particle count is increased, in 
terms of the velocity error norm 
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and the temperature error norm 
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where v and 9 denote the exact solutions for the velocity and temperature. 

The second example is the bar impact problem of Kolsky [33]. The simulations employ a 
linear isothermal equation of state 


P=K (/>/ A> - 1) (83) 

and the parameters shown in Table 3. This problem models the one dimensional motion 
of an elastic bar of length L subjected to a step tensile pressure loading of magnitude P ex t, 
applied at the end x = L at time t = 0. Figures 5 and 6 plot the simulation results for the 
bar midpoint velocity versus time, to a stop time of 0.001 seconds, for two different particle 
counts. The numerical results show good agreement with the exact solution. The results 
shown in Figure 5 are quite similar to those reported by Lu et a. [34], at the same particle 
count, employing an element-free Galerkin method with an explicit integration algorithm. 

The third example problem models the oblique impact of a tungsten alloy (DX2HCMF) 
rod on a steel (SIS 2541) plate, an experiment described in references [29] and [30]. The 


simulations employed a Mie-Gruneisen equation of state [32] and the material properties 
[30,32] listed in Table 4. The cylindrical projectile has a diameter of 0.5 cm and a length 
of 7.5 cm (L/D = 15). The simulation models a 1.5 krn/s impact on a 0.5 cm thick plate 
at a sixty degree obliquity, and was run at three different particle counts. Figures 7 and 8 
show the initial configuration and the simulation results at 100 microseconds after impact, 
while Figures 9 and 10 show sectioned views at 20 and 40 microseconds after impact. This 
simulation illustrates the fact that the present method retains all material fragments and 
models contact-impact of all intact and fragmented material. Table 5 provides simulation 
results, at several particle counts, for the residual rod length and residual rod velocity, 
showing good agreement with the corresponding experimental values (6.38 cm and 1.46 
km/s). The simulation results shown here differ by less than three percent from those 
reported by Lee and Yoo [30] for Lagrangian finite element methods. 

The fourth example problem models the oblique impact of an aluminum sphere on a re- 
inforced carbon-carbon plate. This example does not model a specific experiment; rather it 
was used to measure the improved computational efficiency of the method described here, 
in a hypervelocity impact application of current research interest [6,36]. The projectile ma- 
terial, diameter (0.618 cm), and impact velocity (7 km/s) represent a typical orbital debris 
impact threat [4], while the target plate thickness (0.47 cm) is characteristic of reinforced 
carbon-carbon components of the Space Shuttle thermal protection system [6]. Experimen- 
tal studies of the mechanical properties of reinforced carbon-carbon are in progress; the 
material properties used in the simulations are listed in Table 6. A total of six simulations 
were performed, at three different particle counts, to measure the the computational cost 
of density calculations made using interpolations kernels, as compared to the nonholonomic 
formulation developed in the present paper. Each simulation was run to a stop time of ten 
microseconds, sufficient to perforate the target and include the shock loading and fragmen- 
tation processes of central interest in hypervelocity impact applications. Figures 10 through 
15 show representative simulation results. Table 7 lists the wall clock times and processor 
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counts for the simulations, run on a Linux cluster composed of dual processor (3 GHz) nodes. 
The results show the relatively high computational cost of kernel based density calculations, 
in the present hybrid particle finite element context. Wall clock times are increased for all 
three problem sizes and associated processor counts, by an average factor of one third. This 
result is not surprising, since the finite element related portion of the computation is rela- 
tively inexpensive [22] , while most of the particle related computations are performed in two 
sequential routines: one loops over all neighbor particles to determine the density, the second 
loops again over all neighbor particles to compute the particle interaction forces. The effect 
of introducing the nonholonomic density calculation developed here is to eliminate the first 
of these two routines. Considering the very high computational costs of three dimensional 
shock physics problems, the measured reduction in wall clock time is significant. 

CONCLUSION 

The present paper has formulated a new kernel free particle-finite element method and 
demonstrated its application in three dimensional hypervelocity impact simulations. Unlike 
alternative methods, the formulation is derived without reference to any weighting or in- 
terpolation functions for either the density or the rate of dilatation. The improved method 
introduces a new formulation of the thermomechanical Lagrange equations, one which em- 
ploys internal energies as generalized coordinates. This avoids the requirement to perform 
certain Legendre transforms, in order to express the dependence of the pressure and tem- 
perature on entropy, and hence allows for the use of equations of state in standard form. As 
compared to the previous formulation, the revised method is both simplified and computa- 
tionally more efficient. Applications work in progress is focused on the simulation of orbital 
debris impact effects on spacecraft thermal protection materials [6]. 
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Table 1. Simulation parameters for the wall shock problem 


Ratio of specific heats ( 7 ) 

5 

3 

Specific heat (c„) 

1.0 

Reference internal energy (e D ) 

1.0 

Reference density (p Q ) 

1.0 

Reference temperature (0 o ) 

1.0 

Numerical viscosity coefficient (c 0 ) 

1.0 

Numerical viscosity coefficient (ci) 

0.0 

Numerical conduction coefficient (k 0 ) 

0.1 


Table 2. Error norms for the wall shock problem 


Number of particles 

Velocity error norm 

Temperature error norm 

100 

0.17116 

0.08594 

200 

0.12141 

0.06075 

400 

0.08635 

0.04325 

800 

0.06081 

0.03032 


Tabic 3. Simulation parameters for the bar impact problem 


Length of the bar (A, in) 

100 

Bulk modulus (A", psi) 

30.0 x 10 6 

Applied end loading ( P ext , psi) 

50.0 x 10 3 

Reference density (p D , pci) 

0.73 x 10 “ 3 

Numerical viscosity coefficient (c G ) 

1.0 

Numerical viscosity coefficient (ci) 

0.0 
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Table 4. Material properties used in the long rod impact simulations 


Material property 

Projectile 

Target 

Reference density (g/cc) 

17.6 

7.87 

Shear modulus (Mbar) 

1.45 

0.801 

Reference yield stress (Mbar) 

0.0075 

0.0105 

Strain hardening coefficient 

1.15 

0.177 

Strain hardening exponent 

0.49 

0.12 

Thermal softening coefficient 

1.0 

1.0 

Melt temperature (deg K) 

1,700 

1,723 

Specific heat (Mbar-cm 3 per g-deg K) 

0.143e-5 

0.448e-5 

Plastic failure strain 

1.0 

1.0 


Table 5. Simulation results for the long rod impact problem 


Number of particles 

Residual length 
(cm) 

Residual velocity 
(km/s) 

35,992 

6.08 

1.42 

92,498 

6.47 

1.42 

187,826 

6.59 

1.42 


Table 6. Material properties used in the plate impact simulations 


Material property 

Projectile 

Target 

Reference density (g/cc) 

2.70 

1.58 

Shear modulus (Mbar) 

0.271 

0.0718 

Reference yield stress (Mbar) 

0.0029 

0.000771 

Strain hardening coefficient 

125.0 

2.0 

Strain hardening exponent 

0.10 

1.0 

Thermal softening coefficient 

0.567 

-1.0 

Melt temperature (deg K) 

1,220 

3,840 

Specific heat (Mbar-cm 3 per g-deg K) 

0.884e-5 

0.712e-5 

Plastic failure strain 

1.0 

0.5 
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Table 7. Relative computational costs for the plate impact simulations 


Particles 

Density calculation 

Processors 

Wall clock hours 

Relative cost 

21,334 

nonholonomic 

4 

1.2611 

1.000 

21,334 

kernel 

4 

1.7581 

1.394 

63,253 

nonholonomic 

6 

4.5047 

1.000 

63,253 

kernel 

6 

5.8078 

1.289 

140,070 

nonholonomic 

8 

28.8933 

1.000 

140,070 

kernel 

8 

38.1411 

1.320 
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Figure 1: Wall shock problem, velocity versus position at t = 0.3. 
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Figure 2: Wall shock problem, density versus position at t = 0.3. 
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Figure 3: Wall shock problem, pressure versus position at t = 0.3. 
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midpoint velocity versus time 



Figure 5: Bar impact problem, midpoint velocity versus time for 51 particles. 
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Figure 6: Bar impact problem, midpoint velocity versus time for 101 particles. 



103 



Figure 7: Long rod impact problem, element plot of the initial configuration. 
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Figure 9: Long rod impact problem, sectioned particle-element plot at 20 microseconds after 
impact, color on temperature. 


106 



Figure 10: Long rod impact problem, sectioned particle-element plot at 40 microseconds 
after impact, color on temperature. 
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Figure 11: Plate impact problem, element plot of the initial configuration. 
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Figure 12: Plate impact problem, particle-element plot at 10 microseconds after impact. 
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Figure 13: Plate impact problem, element plot at 10 microseconds after impact. 
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Figure 14: Plate impact problem, sectioned particle-element plot at 10 microseconds after 
impact, color on temperature. 
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Figure 15: Plate impact problem, sectioned element plot at 10 microseconds after impact, 
color on effective plastic strain. 
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